Simulating the behavior of systems that include elastic materials undergoing relatively large deformations (e.g., due to a collision with a rigid body) often involves complex matrix algebra or iterative schemes that often do not efficiently provide an accurate solution within a time frame that would allow for interactive simulations.
For example, the Newton method, and various variations, are often used for simulating nonlinear elastic bodies, such as nonlinear multigrid cycles, projective and position-based dynamics, and shape matching. In a typical Newton scheme, once a linear approximation to the governing equations is computed, the solution to the resulting system is determined using either a direct method or a technique selected from a spectrum of iterative methods.
Direct solvers are typically a fairly straightforward way to solve the system that results from the linearization of the governing equations using a Newton scheme, and can be practical for relatively small problems when direct algebra is not very computational and/or memory intensive. Such direct solvers are often relatively resilient to the conditioning of the underlying problem. Additionally, techniques for directly solving such systems of equations for even relatively large models are available, such as high quality parallel implementations using the Intel MKL PARDISO library. However, such direct solvers typically exhibit superlinear increases in computational complexity as the complexity of the model and/or deformation increases. Even with the benefit of parallelism, such direct solving techniques typically take more time than several iterative schemes, especially if relatively few number of iterations are required to be performed to approximate the solution. Direct solving are also typically limited by the time it takes to access memory to store and retrieve the data used. For example, at the core of many direct solvers are forward and backward substitution routines that carry out a very small number of arithmetic operations for each memory access that is required. This often results in grossly memory-bound execution profiles on modern hardware in which the limiting factor is the time it takes to read to, and write from, memory rather than the time it takes to perform the calculations. This is exacerbated for large models that cannot be fit in cache memory of the processor. Further, each iteration of the Newton method is inherently inexact, providing only a step towards a converged solution. That is, using direct solving techniques often results in perfectly solving an inaccurate linearized approximation of the ultimate solution that is sought.
By contrast, using iterative solving techniques, an approximate solution to the linearized problem can be sought with the understanding that with each Newton iteration the problem itself will often change. For example, such iterative solving techniques can include Krylov methods such as Conjugate Gradient (CG), Multigrid, and fixed-point iterations such as Jacobi, Gauss-Seidel (GS), and successive over-relaxation (SOR). The primary benefit of such iterative techniques is that often each individual iteration can be performed relatively quickly, which can allow a user the option to either iterate as much as they can afford in a given time, or alternatively truncate the iterative process when the approximate solution is acceptable (e.g., each iteration is producing a solution that changes less than a threshold amount from the previous iteration). Additionally, many iterative solving techniques are assembly-free, alleviating the need to construct or store a stiffness matrix. In fact, some of the most efficient techniques go to great lengths to minimize memory footprint while leveraging single instruction, multiple data (SIMD) and multithreading techniques.
However, iterative solving techniques often have different challenges. For example, local techniques such as Jacobi, GS, and SOR are slow to capture global effects, as they propagate information at a limited speed across the mesh of a model. As another example, many Krylov methods prioritize the most important modes that contribute to a high residual. In a more particular example, considering a system with a few tangled elements that create large local forces, elements suffering from small errors will be relatively neglected by a technique such as Conjugate Gradients, while the solver focuses computational resources on the highly tangled elements before resolving the bigger picture. While multigrid often has performance advantages, it can be relatively difficult to configure to work robustly, and might be less appropriate for certain types of objects such as thin elastic objects (e.g., a thin flesh layer on a simulated face). Preconditioning can accelerate the convergence of iterative solvers, but—in contrast to certain fluids simulation scenarios—the accelerated convergence may not justify the increased per-iteration cost. Preconditioners based on incomplete factorizations are typically memory bound as they require matrix assembly, and generally require an expensive re-factorization at each Newton iteration. Further, the same factorization overhead is incurred in each Newton iteration regardless of how closely the Newton method is to convergence. For example, in a case where the Newton method is nearly converged such that just a handful of iterations would suffice to solve the linearized equations, the same factorization overhead is incurred which can contribute a relatively large portion of the total time to solve the linearized equations. Multigrid-based preconditioners can achieve more competitive performance in some circumstances, and have been primarily explored in the area of fluid simulation rather than in the simulation of nonlinear deformable solids.
When fidelity and realism of the simulation is a goal, physics-based methods are often employed, such as the Finite Element Method, which has been used to animate a diverse spectrum of behaviors. Grid-based, embedded elastic models have also been used frequently due to their potential for performance optimizations, and can also be used with shape-matching approaches. Such grid-based models form the foundation for a class of relatively efficient, multigrid-based numerical solution techniques. Various techniques to accelerate simulation performance have been proposed. For example, using optimized direct solvers, delayed updates to factorization approaches, leveraging the Boundary Element Method to approach real-time deformation, and similar formulations that abstract away interior degrees of freedom to accelerate collision processing. However, while the techniques described above may improve performance, they may suffer from significant disadvantages, such as being limited by memory access, superlinear increases in complexity with model size, delays in capturing global effects, etc.
Accordingly, new systems, methods and, media for simulating deformations of nonlinear elastic bodies are desirable.